Plasticity and Dislocation Dynamics in a Phase Field Crystal Model 
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The critical dynamics of dislocation avalanches in plastic flow is examined using a phase field 
crystal (PFC) model. In the model, dislocations are naturally created, without any ad hoc creation 
rules, by applying a shearing force to the perfectly periodic ground state. These dislocations diffuse, 
interact and annihilate with one another, forming avalanche events. By data collapsing the event 
energy probability density function for different shearing rates, a connection to interface depinning 
dynamics is confirmed. The relevant critical exponents agree with mean field theory predictions. 
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Materials yield and deform plastically under large ex- 
ternal stress. While the yield surface and the plastic flow 
are well described by various continuum theories [TH3], 
what happens microscopically during a plastic deforma- 
tion is still not fully understood. On atomic scales, ex- 
ternal stress is carried by localized crystal defects, such 
as dislocations and disclinations. They are created un- 
der stress and interact with each other. Although the 
properties of individual defects and their interactions are 
well known [H [5], their collective behavior under exter- 
nal stress is complicated. It gives rise to scale invariant, 
power-law distributed phenomena [BUT?], in strong re- 
semblance to the scaling behavior near a critical point. 
These phenomena include dislocation slip avalanches of a 
broad range of sizes, and acoustic emission with a power 
law power spectrum. 

Recently, evidence has accumulated that these scaling 
phenomena reflect an underlying non-equilibrium critical 
point [6] [7], i.e. a point where the non-equilibrium steady 
state fluctuations are governed by a diverging correla- 
tion length |13j . as appears to be the case in magnetic 
materials pH] and even turbulence [TS]. Such a critical 
point would exist at the boundary between two distinct 
regimes, one of which would be a glassy, activated regime, 
the other would be a genuine plastic flow regime. Up to 
now, experimental and computational data have focused 
on the glassy regime: for example, Weiss et al. have 
measured the acoustic emission signal from creep defor- 
mation experiments on single crystal ice and found that 
the event size distribution follows a power law over 4 
decades [3 [8]. Miguel et afs dislocation dynamics simu- 
lations in two dimensions showed event size distributions 
following a power law, with a rate-dependent cutoff, over 
approximately 2 decades [5] . Zaiser et al. have reported a 
data collapse of the event size distribution with different 
external stresses [Sj. 

The purpose of this Letter is to approach plastic defor- 
mation from the flow side of the non-equilibrium critical 
point, manifested in the strain-rate dependence of the 
acoustic emission. Importantly, we are able to system- 
atically vary the strain-rate in simulations, and more- 
over we relate the critical point underlying plastic flow 
at finite strain rates to the scaling of magnetic domain 



wall depinning [6] [HI [17]. We find remarkable agree- 
ment between simulations and analytical mean field the- 
ory predictions of exponents, and in addition arc able 
to show that the strain-rate data exhibit collapse. Our 
results strongly support the critical point picture of plas- 
ticity, and suggest new experiments. We study disloca- 
tion avalanches during plastic flows using the phase field 
crystal (PFC) model[18j[TH]. This approach is well-suited 
to this problem, because it can be performed at finite 
temperature, for large systems, and over long time pe- 
riods. The PFC model describes the dynamics of the 
local crystalline density field, and has been shown to 
give an excellent account of numerous materials proper- 
ties including polycrystalline solidification, vacancy dif- 
fusion, grain growth, grain boundary energetics, epitaxial 
growth, fracture 19J , grain coarsening[20 , elasticity 21 , 
dislocation annihilation, glides and climb 22J, as well as 
vacancy dynamics [23] < The model can be derived from 
density functional theory and extended to the case of bi- 
nary alloys 24 . In this paper, we augment the model 
to treat external shearing forces by adding an advective 
term to the dynamics near the boundary. By adjusting 
the shearing force and measuring the resulting avalanche 
statistics, a data collapse is obtained, that is consistent 
with proximity to a domain wall depinning point at finite 
temperature [rJlUT]. 

The Model:- The phase field crystal (PFC) model is 
given by the free energy density [HI [19] 



f=^ 2 + l?P+- 2 P 2 + , 
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where r is the undercooling and p(x, t) is the local den- 
sity. The dynamics associated with this free energy is 
conservative, relaxational and diffusive, and systemati- 
cally derivable from density-functional theory [25]. In or- 
der to study the plastic response of the PFC model under 
shear, we add a shearing term to the dynamical equation: 
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is the shearing profile and vo is the magnitude of the 
shearing, A is the penetration depth, a and /3 control 
the range and time scale respectively of elastic inter- 
actions (phonon excitations) propagating through the 
medium [2 1, F = J f(x)d d x is the total free energy and r\ 
is the thermal noise satisfying the fluctuation-dissipation 
theorem (r)(x,t)r]{x",t')) = -e\7 2 5{x - x')5(t - t'). Here 
e is the noise amplitude. It is directly proportional to 
the temperature, ksT. The value of vo controls the mag- 
nitude of the shearing force; the penetration depth, A, 
controls how deep the shearing force extends into the 
material. In all simulations we set A <C L yi so the ac- 
tual value of A does not affect our simulation results. 
The boundary conditions are periodic in x and fixed at 
y = 0, L. The PFC model we used allows propagating 
sound modesJST] when /3 = 0.9; as long as (3 is non-zero 
and 0(1), we do not expect that our scaling results will 
be sensitive to its precise value. Values of that are 
O(10) would correspond to very over-damped dynamics, 
and could model viscoelastic behavior [2"T] that is outside 
the scope of our work. 

One of the advantages of using the PFC model is that 
we do not have to impose any ad hoc assumptions about 
the creation and annihilation of dislocations. Recall 
that in dislocation dynamics simulations, dislocations are 
treated as elementary particles and usually only the far 
field interaction between dislocations is captured. When 
dislocations get too close to each other (a few atomic 
spacings) , the highly nonlinear interaction between them 
is not captured and more importantly, the annihilation of 
dislocations is not accounted for. The standard practice 
is then to impose some annihilation rules — dislocations of 
opposite topological charges annihilate when they get too 
close to each other[S]. Similarly, dislocations have to be 
created by hand when the local strain is high. Although 
these rules are consistent with our physical intuition, par- 
ticular ways of implementing them are sometimes diffi- 
cult to justify. However, because the PFC model captures 
the nonlinear elastic behavior of a crystal, the interaction 
between dislocations is completely captured. In addi- 
tion, because the PFC model simulates the atoms in the 
lattice (the PFC density is periodic in its ground state 
with peaks representing atoms and troughs inter-atomic 
space), but not the dislocations themselves, creation and 
annihilation of them are also naturally captured as col- 
lective excitations of the lattice. No ad hoc rules or as- 
sumptions have to be imposed. 

We solved Eq. Q in a 2D rectangular domain. The 
crystal under shear is initially perfectly triangular. As 
the crystal is sheared, dislocations are created near the 
fixed boundaries y — 0, L where the stress is higher. They 
then propagate into the bulk. They interact with each 
other and form avalanches. To quantify the avalanche 
activity, we calculate the total speed of all dislocations 
in the domain, V(t) — J2i=i K»l> where Ndi S (t) is the 
number of dislocations in the system at time t and u; is 
the velocity of the i-th dislocation. This measure is simi- 
lar to the acoustic emission signal in Weiss et al.'s single 



crystal ice experiments. As dislocations are generated 
and interact with each other in the domain, in addition 
to the fast avalanching dynamics, quasi-static structures, 
such as grain boundaries, can form. These slow dynamics 
should not be measured because they are really not part 
of the avalanches. This leads to the distinction between 
fast-moving and slowly-moving dislocations introduced 
by Miguel et al. 9 . In essence, they introduced a cutoff 
in dislocation speed and measure only dislocations with 
speed higher than the cutoff. In that way, they tried to 
retain only the avalanche activities in the acoustic emis- 
sion signals. 

We employed a different method to tackle this problem. 
Instead of simulating a very large system, with all sorts 
of dislocation activities, we simulated a moderate size of 
system with approximately 10000 atoms. For this system 
size, dislocation avalanches come and go, i.e., not many 
dislocations are left in the system after every avalanche. 
As a result, no grain boundaries, or slow dynamics, is 
present and we obtain clean avalanche data. It is fair to 
mention that this method severely limits the system size, 
and thus the resulting avalanche sizes. The system size 
we chose contains approximately 100 dislocations in the 
largest avalanche events. The tradeoff, which we exploit, 
is the cleanness of the avalanche signal and the speed 
of the resulting simulations. Different methods, such as 
those we mentioned above, would have to be employed if 
larger avalanche sizes are desired. 

We count the number of nearest neighbors of each 
atom, n.i, using the Delaunay triangulation method in 
computational geometry [2l)I |2"7] . Because we have n; = 6 
for every atom in a perfectly triangular crystal, and be- 
cause there are no vacancies introduced in the version 
of the PFC model simulated here (vacancies can be in- 
troduced into the PFC model by breaking the up-down 
symmetry of the PFC free energy, as detailed in |23j). 
any atom having rii 7^ 6 is sitting next to a dislocation. 
Because the PFC exhibits emergent rigidity in the region 
of the phase diagram studied here, these 'defect atoms' 
essentially track the locations of dislocations. Instead of 
measuring the total sum of dislocation speeds, V(t), we 
then measure the total sum of these defect atoms' speeds, 



N(t) 

v(t) = £ 



(4) 



where N(t) is the number of defect atoms and Vi is the 
velocity of defect atom i. Note that the velocity of a de- 
fect atom is not the velocity of any atom in the system, 
but the velocity of the dislocation it is tracking. Because 
the two measures, V(t) and V(t) are proportional to each 
other with the proportionality constant being the mean 
number of defect atoms sitting next to a dislocation, we 
can use the latter for convenience. Fig. [T] shows the 
typical time dependence of N(t) from a simulation with 
parameters dx — 37r/8, dt = 0.025, L x = L y = 512, 
(a) 2 = 255, /? = 0.9, v = 1.581, p a = 0.3, e = 1.5, 
A = 40.0 and r = —0.5. N(t) changes as dislocations are 
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FIG. 1: The number of dislocations in a sheared PFC crystal. 
Intermittent events with sizes differing in orders of magnitude 
is observed. Parameters are (a) 2 = 255, /3 = 0.9, vo — 1.581, 
po = 0.3, e = 1.5, A = 40.0 and r = -0.5. 



FIG. 3: (Color online) The probability distribution of the 
event energy during dislocation avalanches, for different val- 
ues of the shearing rates. 




FIG. 2: The total speed of defect atoms in a sheared PFC 
crystal. Intermittent events with sizes differing in orders of 
magnitude is observed. Parameters are (a) 2 = 255, /3 = 0.9, 
v = 1.581, po = 0.3, e = 1.5, A = 40.0 and r = -0.5. 



being created and annihilated. There are intermittent 
events of creation of dislocations, with number of dislo- 
cations involved ranging from a few to 80. Fig [2] shows 
the acoustic emission signal, V(t), in the same simula- 
tion. Similar to N(t), the signal ranges from to 400, 
with intermittent pulses of various sizes. In order to mea- 
sure the avalanche event size, we introduce a cutoff to the 
signal V cu t — 15 for shear velocity vq — 0.765. The cutoff 
is increased proportionally to the shear velocity, and is 
primarily used to define the avalanche size, rather than 
to remove slowly moving dislocations from the analysis. 
The signal is then partitioned into individual avalanche 
events. The probability distribution of the event energy, 



E 



V 2 {t)dt, 



(5) 



where tbegin and t en d are the starting and ending time 
of the event respectively, can be measured. For cutoff 
values small compared to the signal (i.e. V cu t = 15, cor- 



responding to the activity of 3-4 dislocations due to ther- 
mal creep) , the result is insensitive to the cutoff. For each 
shearing rate, at least 10 different realizations are run to 
obtain a statistically meaningful result. This results in 
about 8000 avalanche events for each shearing rate. Fig. 
[3] shows the event size distribution for different shearing 
rates. We find that the distribution follows a power law 
for small event sizes and cuts off at larger sizes, with 
the cutoff size depending on the shearing rate. The data 
is somewhat noisier towards the end of large event size 
because large events are rare. 

Scaling Behavior of the avalanches:- Analogous to the 
scaling behavior in models of crackling noise |14j . we 
propose that there is a non-equilibrium critical point, 
vq = v c , in the system and we expect, as E — > 00, the 
data around the critical point to collapse in the form 



P{E,v) ~ E~ T f(Ev^ 



(6) 



where P{E, v) is the probability distribution of event 
energy, E, v = 1 — Vq/v c is the reduced shearing rate 
with vo being the shearing rate and v c being the critical 
shearing rate, r and /j, are two critical exponents. As 
v — > 0, P(E,v) tends to a power law P(E,v) ~ E~ T . 
Fig. [4] shows an attempt to collapse the data, using 
the equivalent scaling form P(E,v) ~ v T ^g{Ev^), with 
t = 1.5, /i = 2 and v c = 1.5 and universal scaling func- 
tion g(x) = x~ T f(x) shown in the collapse. Logarith- 
mic binning is performed and singletons are ignored to 
obtain P(E,v). We obtain a satisfactory data collapse 
over 4 decades, with v ranging from 0.15 to 0.49. As 
E — > 0, the collapse function g{x) approaches the power 
law g{x) ~ x~ 3 ' 2 , which agrees with f(x) — > const in 
Eq. ([6]) for x — > and r = 3/2. The collapse constrains 
the numerical exponent to the range r = 1.5 ± 0.2. This 
agrees with the experimental result of t = 1.5 [55]. Sim- 
ilarly, for adiabatic stress increase in the pinned regime, 
Zaiscr finds r = 1.4 [6]. Dimiduk reports r = 1.5 — 1.6 
in experiments at fixed compression stress that leads to 
shearing |29j . In contrast to the universal distribution 
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FIG. 4: (Color online) Data collapse of the probability dis- 
tribution of the event energy during dislocation avalanches, 
with t = 1.5, /J. = 2 and v c — 1.5. 



of time integrated avalanche signals discussed here (see 
eq.([5])), the signal amplitude statistics likely depend on 
details of the system: The distribution of energy ampli- 
tudes decays with the exponent r = 1.8 [9] in simulations, 
and with r = 1.6 in experiments |5|. The accoustic emis- 
sion intensity exponent is r = 1.8 in simulations at fixed 
stress [3"U] . 

Studies at adiabatically slow shear rate have suggested 



analogies to the depinning transition of magnetic domain 
walls [H [17l |3T] , although one exponent that may deviate 
in simulations is discussed in [32]. The exponents found 
in our collapse are consistent with the domain wall de- 
pinning picture. As argued previously, mean field theory 
(MFT) is expected to give exact scaling results in this 
case. The MFT values for the exponent r is r = 1.5 
O [T7j . The MFT value for the exponent pi can be cal- 
culated from the MFT of QU QH |3J] to be fj, = 2. These 
MFT values for the exponents r and fj, lead to a satisfac- 
tory collapse of the numerical avalanche size distributions 
at different shear rates. At zero temperature the critical 
shear rate is v c — 0. Our simulations, however, are per- 
formed at finite temperature T. Temperature-induced 
dislocation creep causes the critical shear rate to appear 
to be nonzero, with the apparent v c (T) — > as T — >• 0. It 
also causes the scaling collapse in Fig. [4] to be less precise 
than in zero temperature studies. We hope to report on a 
theoretical investigation of the temperature dependence 
and a comparison with experiments at finite shear rate 
in a future publication. 
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